Stochastic dynamics for a single vibrational mode in molecular junctions 



O 

(N 

C 



A. Nocera 1 , C.A. Perroni 2 , V. Marigliano Ramaglia 2 and V. Cataudella 2 

1 Dipartimento di Fisica E. Amaldi, Universita' di Roma Tre, 

Via della Vasca Navale 84, 1-0014-6 Roma, Italy 

2 CNR-SPIN and Universita' degli Studi di Napoli Federico II 
Complesso Universitario Monte SantAngelo, Via Cintia, 1-80126 Napoli, Italy. 

We propose a very accurate computational scheme for the dynamics of a classical oscillator coupled 
to a molecular junction driven by a finite bias, including the finite mass effect. We focus on two min- 
imal models for the molecular junction: Anderson-Holstein (AH) and two-site Su-Schrieffer-Heeger 
(SSH) models. As concerns the oscillator dynamics, we are able to recover a Langevin equation 
confirming what found by other authors with different approaches and assessing that quantum ef- 
fects come from the electronic subsystem only. Solving numerically the stochastic equation, we 
study the position and velocity distribution probabilities of the oscillator and the electronic trans- 
port properties at arbitrary values of electron-oscillator interaction, gate and bias voltages. The 
range of validity of the adiabatic approximation is established in a systematic way by analyzing the 
behaviour of the kinetic energy of the oscillator. Due to the dynamical fluctuations, at intermediate 
bias voltages, the velocity distributions deviate from a gaussian shape and the average kinetic energy 
shows a non monotonic behaviour. In this same regime of parameters, the dynamical effects favour 
the conduction far from electronic resonances where small currents are observed in the infinite mass 
approximation. These effects are enhanced in the two-site SSH model due to the presence of the 
intermolecular hopping t. Remarkably, for sufficiently large hopping with respect to tunneling on 
the molecule, small interaction strengths and at intermediate bias (non gaussian regime), we point 
out a correspondence between the minima of the kinetic energy and the maxima of the dynamical 
conductance. 
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I. INTRODUCTION 



In recent years it has become possible to fabricate 
electronic devices where the effective element of a junc- 
tion is a single molecule placed between two metallic (or 
semiconductor) electrodes* 1 - - — Due to the small size of 
molecules, the charging of the molecular bridge is often 
accompanied by significant changes of the nuclear geom- 
etry, indicating a strong coupling between electronic and 
nuclear (in particular vibrational) degrees of freedom. 
For example, some authors^ have recently proposed a the- 
oretical explanation of the switching mechanism, actually 
observed in different molecular junctions,^, based on the 
electron-phonon interaction. Understanding and control- 
ling the effect of this interaction onto the electric current 
through molecular devices is not only important for the 
field of molecular electronics but establishes a strong link 
also to the physics of Nano-ElectroMechanical Systems 
(NEMS)i^ii Recent experiments show the possibility to 
use single electron transistors coupled to a mechanical os- 
cillator as high sensitive position^ and mas*^ sensors. 
One of the challenges is to understand, control and use 
the interplay between a quantum v v detector ' ' (electron 
transistor) and a classical mechanical system. 

The simplest molecular conduction junction comprises 
two metallic electrodes connected by a single molecule. 
Such a junction, including the effect of electron-phonon 
interaction, can be described by the Anderson-Holstein 
(AH) modelfii The molecule is represented by one elec- 
tronic level interacting linearly with a local vibrational 
degree of freedom and connected through tunneling with 
free-electron metals. Electron transport within this 



model has received a lot of theoretical attention! 12 ' 13 De- 
spite the conceptual simplicity, it gives rise to a very rich 
physics. Several approaches have been adopted depend- 
ing on the relative energy scales in the problem. When 
the characteristic frequency of the oscillator ujq is of the 
same order of magnitude of the tunneling frequency of 
the electrons on the molecule (~ T), a quantum treat- 
ment of oscillator dynamics is necessary. In this case, it 
is useful to consider separately the limits of weak and 
strong electron-phonon interaction strength relative to 
the coupling of the level to the leads. The former cor- 
responds to nonresonant phonon-assisted electron tun- 
neling, mostly encountered in experiments in inelastic 
electron tunneling spectroscopy 14 (IETS), and theoret- 
ically understood within Non Equilibrium Green Func- 
tion (NEGF) formalism*^ - — In the case of stronger ef- 
fective electron-phonon coupling, the perturbative treat- 
ment breaks down, the conduction shows phononic block- 
ade at small bias (Franck-Condon effect 1 ^—) and one can 
observe phonon sidebands in the conductance spectra^ 2 - 
The interpolation from weak to strong electron-phonon 
coupling regime in a full quantum description of electron 
and phonon subsystems is a very challenging problem. 
However, when the vibrational motion is slow compared 
to the electronic tunneling rate (uq « T), one can apply 
a simplified scheme in the spirit of Born-Oppenheimer 
approximation. Indeed, in the zero-order static theory 
(adiabatic ratio ujq/T t— > 0, oscillator mass m i— > oo), one 
neglects the kinetic energy of the oscillator and obtains 
an exact electronic problem where the oscillator position 
enters as a parameter to be determined self-consistently. 4 
Some authors 2 ^— have already considered the possibil- 
ity to construct corrections to this picture, within the AH 
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model, in different parameter regimes and with different 
techniques. Using the action functional formalism on the 
Keldysh contour for the full interacting electron-phonon 
problem, Mozyrsky et al. 2 ^ obtain an effective Langevin 
equation for the oscillator field in the limit where the 
electron-leads are considered as zero temperature ther- 
mostats. It comprises a position-dependent dissipation 
term and white noise force. In the strong electron- 
oscillator coupling regime, where the model shows bista- 
bility, they find that the oscillator field acquires an ef- 
fective temperature linearly related to the bias, Vuasi 
as a consequence of the coupling to the electronic bath. 
Pistolcsi et al. 2A starting from the Langevin equation, 
generalize the previous results solving numerically the 
Fokker-Plank equation associated with it, focusing again 
on the extremely strong coupling regime. They obtain 
the dependence of the current on the transport and gate 
voltages, as well as address the problem of mechanical 
switching between the metastable states of the oscillator 
potential. 

A similar approach, based on the Feynman- Vernon in- 
fluence functional, was recently adopted by Hussein et 
al2^ obtaining the same Langevin equation. However, 
they do not solve this equation and study the phase space 
portraits of the Newton equation obtained in absence of 
electronic bath induced noise. Moreover, they extend 
this analysis to a molecular double dot molecular Hamil- 
tonian. For both cases they have explained the features 
of effective potential and friction terms entering the os- 
cillator equation of motion. 

The approaches discussed above, based on an expan- 
sion of the action in the adiabatic limit, are not able 
to disentangle the origin of the quantum effects in the 
Langevin equation for the oscillator. We propose here an 
alternative and more direct way that allows us to clarify 
this point. 

We construct systematically the dynamical finite mass 
corrections to the static theory and their influence on 
the transport problem in the following way: we per- 
form an adiabatic expansion on the electron-oscillator 
self-energy following the approach of Ref.26, obtaining 
a corresponding expansion of the Green function. This 
expansion gives rise to the same friction and fluctuat- 
ing terms obtained with action functional techniques, but 
clarifies that the quantum effects v hidden ' in the stochas- 
tic equation come only from the electronic subsystem. 

We numerically solve the Langevin equation deriving 
the position and velocity distribution probabilities of the 
oscillator for a very large range of the relevant parame- 
ters. We find that, at intermediate bias voltages, the ve- 
locity distributions P{v) deviate from a gaussian shape 
as a result of the coupling of the oscillator with the out- 
of-equilibrium electronic bathj^l Correspondingly, the ki- 
netic energy of the oscillator shows a non monotonic be- 
haviour as function of the bias due to the slight change 
of the force exerted on the oscillator. 

We study transport properties like the current-voltage 
characteristic and the conductance, observing in the AH 



model a dynamical reduction of the the * polaronic ' shift 
and the broadening of the electronic resonance due to 
the average over the nonequilibrium position distribu- 
tion probability of the oscillator. We note an interesting 
strong enhancement of current in the non gaussian region 
at intermediate bias, where the infinite mass approxima- 
tion prescribes very small conduction. 

It is of paramount importance to study the range of 
validity of the adiabatic approach. Making a thorough 
investigation of this issue is crucial in order to get a 
match with experimental results and exact theoretical 
calculations. We establish the range of validity of the 
adiabatic approximation analyzing systematically the be- 
haviour of the oscillator's kinetic energy. We are able to 
build up a diagram for the validity of classical approxi- 
mation, identifying Quantum (QR), Classical- Adiabatic 
(CAR) and Classical Non-Adiabatic (CNAR) Regions. 
We compare the classical kinetic energy of the oscillator 
with the Debye temperature (fc^Tc ~ Tiujq) to distin- 
guish between QR against CAR regimes, and with elec- 
tron energy scale (~ TiT) to distinguish between CAR 
against CNAR regimes. 

We extend this analysis to the case of a molecular 
Hamiltonian composed by a couple of sites interacting 
with a single vibrational mode in the SSH modeler— In 
this case, because of the direct coupling of the electron- 
oscillator interaction to the intermolecular hopping, one 
expects that the role of the dynamical fluctuations be- 
comes crucial to determine the correct features of the 
observables inherent to the transport problem. 

In the limit of symmetric coupling with leads, we are 
able to construct again a Langevin equation for the oscil- 
lator dynamics, very similar to that derived in AH model. 
In this case, it is possible to study the effect of the dy- 
namics of the classical vibrational mode on the electron 
hopping through the two molecular sites. Vice- versa, one 
can also study the effect of intermolecular electron degree 
of coherence onto the vibrational dynamics. 

The new intermolecular electronic hopping scale t in- 
troduces a reduction of the CAR in the validity diagram. 
This reduction becomes important if t/hT > 1 with the 
occurrence of QRs. These new features are due to the 
stronger nonmonotonic behaviour of the kinetic energy 
of the oscillator as function of the bias voltage. For suf- 
ficiently large t/hT, small interaction strengths and at 
intermediate bias (non gaussian regime), the kinetic en- 
ergy curves show well defined minima. We point out a 
correspondence between these minima and the maxima 
of the conductance. Again, as concerns the electronic 
transport properties, the dynamical fluctuations favour 
the conduction far from the two electronic (static) reso- 
nances. In the SSH model, as already stressed above, the 
effects of the dynamical fluctuations become even more 
important. One can observe again a complete erasing of 
the bistability and hysteretic behaviour predicted by the 
infinite mass approximation. 

The paper is organized as follows: In Sec. II we present 
the single level case within the AH model. In Sec. Ill we 
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deal with the case of two-site-SSH model. 



II. THE ANDERSON-HOLSTEIN (AH) MODEL 

The spinless Anderson-Holstein model is the simplest 
model of a molecular junction including the effect of 
electron-phonon interaction. The molecule is modeled 
as a single electronic level interacting locally with a 
single vibrational mode. The electronic system is de- 
scribed by the standard junction Hamiltonian T-L e i = 

H mo l + Hfun + Hi ea ds , with 

H mol = E g Sd, (1) 



(r*:,a | — >■ Tq,). In this paper we will measure length in 
units of Xq = and energy in units of UT. Finally, the 

leads will be considered as zero temperature thermostats. 

In the next subsections, we will first (subsection A) an- 
alyze the coupled electron-oscillator problem in the limit 
of infinite mass for the oscillator. We will then indicate 
how to construct (subsection B) the stochastic Langevin 
equation for the dynamics of the oscillator including the 
finite mass effect. In the subsection C we will solve nu- 
merically the stochastic equation and analyze the effects 
of the oscillator dynamics on the electronic observables 
inherent to the transport problem (I-V characteristic and 
conductance) . 
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The molecular electronic level has energy E g and 
creation (annihilation) operators ^(d). The operators 
c\ a (ck, a ) create (annihilate) electrons with momentum 
k and energy 6k, a — £,k,a — Ha in the left (a — L) or right 
(a = R) free metallic leads. The chemical potentials in 
the leads, fiL and fiR> are assumed to be biased by an 
external voltage, eVuas = Ml — HR- Electronic tunnel- 
ing between the molecular dot and a state in the lead 
has amplitude Vk, a - We consider the oscillator dynam- 
ics " classical ' from the beginning and described by the 
position and momentum variables x,p. 

The Hamiltonian of the oscillator is given by 



A. Out of equilibrium Born-Oppenheimer 
approximation: infinite mass (static) case 

When the vibrational motion of the oscillator is slow 
with respect to all electronic energy scales, it is possible 
to decouple oscillator and electronic dynamics. In the 
spirit of Born-Oppenheimer approximation, we consider 
the limit m \— > oo in the full Hamiltonian disregarding 
the kinetic energy of the oscillator. The electronic dy- 
namics is therefore equivalent to a non-interacting res- 
onant single level problem with energy level renormal- 
ized by the 'polaronic' shift E g M- E g + Xx. The re- 
tarded (advanced) Green functions G r w (uj, x) and the 
lesser (greater) Green functions (uj, x) in stationary 

nonequilibrium conditions are derived within the Keldysh 
formalism 16 ' 17 and depend parametrically by the dis- 
placement coordinate x. Starting from the force exerted 
on the oscillator 



F = -mugx + \(n e i)(x), 



(6) 



H SC 



2m 



(4) 



characterized by the frequency luq and the effective mass 
m. The interaction (typically of electrostatic origin) is 
provided by a simple linear coupling between the electron 
occupation on the molecule, cftd, and the displacement of 
the oscillator 



H,, 



Xxd) d, 



(5) 



where A is the electron-oscillator coupling (EOC) 
strength. The overall Hamiltonian is therefore given by 

it = T~Lel + H osc + Hint- 

In the following, the coupling between the electron sys- 
tem and the vibrational mode will be often described 
in terms of the electron-phonon coupling energy E p — 
A 2 /(2mcjQ), while the coupling to the leads by the tun- 
neling rate i\ Q = 27r / 9 Q |Vfc j „| 2 /7i (the full hybridization 
width of the molecular orbital is then hTk = K^k,L + 
K^k,R), where p a is the density of states of the lead a. 
For the sake of simplicity, we will suppose flat density of 
states for the leads within the wide-band approximation 



where 



(n el )(x)= J2 ^ 2 r« /I^HIG^)! 2 , (7) 

T r> J 



a=LM. 

with / a (w) Fermi function of the lead a = R,L, one 
can therefore straightforwardly compute the expression 
of the generalized potential in nonequilibrium condi- 
tions (obtained applying symmetric bias unbalance Hr = 

-eVbias/2, Ml = eVbi as /2) 



U(x) 



arctan 



f 



2„2 



-mujr.x 



Xx 
T 



E 

a=L,R 



Eg — XX 



2tt 



Eg , Xx ) - — ln[4( Ma - E q - Xx) 2 

nr/2 J 8tt l v ^ 9 ' 



-{nvf 



(8) 



This generalized oscillator potential depends paramet- 
rically by the spring constant maij = k, the EOC 
strength A, the energy of the electron level E g (which 
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FIG. 1. (Color online). Spatial dependence of the dimension- 
less generalized static potential U(x) (panels A,B), friction 
coefficient A(x) (panels C,D), fluctuating term D(x) (pan- 
els E,F) for symmetric E g ~ E p and asymmetric E g < E v 
minima and different values of bias, Vuaa = (solid (black) 
curve), Vbias = 2 (dashed (red) curve), Vuas = 4 (dashed dot 
(green) curve). The potential is expressed in hT units (U — 
U/HT), the friction coefficient in muo units (A = A/muo), 
the fluctuating term in A 2 /cjo units, (D — D/{\ 2 /wo))- Vbias 
values are expressed in fiT/e units, where e is the electron 
charge. The dimensionless position variable x is defined as 
x — x/xa with xa = — 

' mi.)- 



can be considered a gate potential), the coupling to the 
leads T and finally by the bias, Vbias- In FigQ] (panels 
A,B), we present some features of the generalized poten- 
tial U(x) in the strong coupling regime (E p > hT), where 
the potential shows several minima. For E g ~ E p and 
not too large bias (panel A, Fig[T]), the potential devel- 
ops two symmetric minima near x ~ (corresponding to 
(n e i) ~ 0) and x ~ —1 (corresponding to (n e i) ~ 1) sep- 
arated by a barrier whose height is roughly proportional 
to E p . This bistable regime corresponds to the physical 



situation where the bare electron level E g is above the 
chemical potential of both leads, while the renormalized 
charged level E g — 2E p is below them. The molecule 
can stay in one minimum or in the other. If we increase 
the bias Vbias , the potential U(x) shows a third mini- 
mum corresponding to average electron occupation on 
the molecule (h e i) — 1/2 and, for sufficiently large Vbias, 
only this minimum remains. If E g < E p the potential 
also shows two or more minima but they are asymmetric 
(panel B, FigJT]). For sufficiently large bias, the common 
feature is the existence of a single minimum correspond- 
ing to occupation (h e i) ~ 1/2. 

In the above analysis the displacement x has been used 
as a free parameter. Actually, the only x values relevant 



for the electronic properties in the static approximation 
are those which solve the equation F(x) = 0. These solu- 
tions depend parametrically by all the parameters of the 
theory (in particular by the bias Vbias)- This may yield 
transitions between different local minima in the poten- 
tial, determining in the electronic current-voltage char- 
acteristic the onset of interesting non linear phenomena 
like hysteresis, bistability and Negative Differential Re- 
sistance (NDR)4 Indeed, the authors of Ref.4 proposed 
a polaron mechanism within the AH model to explain 
such phenomena, effectively observed in transport ex- 
periments on molecular devices. However, the results 
of the static approximation can be strongly modified by 
dynamical effects. Indeed, corrections due to the finite 
(though large) mass of the oscillator are expected to be 



important 
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As we shall see in the next sections, the 



inclusion of the finite mass effect on the oscillator dy- 
namics gives rise to a stochastic Langevin equation with 
a position dependent dissipation term and white noise 
force. The stochastic fluctuations of the oscillator mo- 
tion will strongly modify the current-voltage characteris- 
tics obtained in the infinite mass approximation. 



B. Dynamical (finite mass) corrections to static 
case: setting Langevin equation for the oscillator 

Within the static approximation (infinite mass), the 
main effect of the nonequilibrium fast electronic environ- 
ment is the modification of the force (Eq.©) experienced 
by the mechanical oscillator. The oscillator has no dy- 
namics at all and the displacements x are " frozen ' in 
suitable points of the configuration space given by the 
equation F(x) = 0. In this section, we show how to 
include dynamical corrections due to finite mass of the 
oscillator. 

First of all, we should include the time dependence of 
the oscillator dynamics in the Hamiltonian of the elec- 
tronic problem (where we still neglect the oscillator ki- 
netic energy, 7i i-> H(x(t))). Using the extension of the 
Keldysh formalism to time dependent cases^ we can 
solve the Dyson and Keldysh equations for the molecu- 
lar Green functions which now depends on times t and t' 
separately. The retarded molecular Green function can 
be obtained analitically (with E g (t) — E g + Xx(t)) 



G r (t,t') 



-9{t-t')e 



-»r/2) 



(9) 



and depends in non linear way by the entire dynamics 
x(t) of the oscillator. Indeed, the dependence by the 
entire history of the oscillator motion prevents us from 
Fourier transforming the Green function and leads to an 
intractable problem. 

In order to overcome this difficulty, we resort to an 
adiabatic approximation of the molecular Green func- 
tion. As thoroughly discussed in the Appendix A, one 
should restart from the Dyson equation for the Green 
function (Eq. (|Aip ) and perform an expansion in the 
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electron-phonon self-energy where a separation between 
" slow ' and v fast ' time scales was preliminary accom- 
plished fEqs. (|A3|A4j) ). The expansion, performed with 
respect to the "slow' time (t + t')/2, allows to disentan- 
gle the non-local time dependence of the Green function. 
In the end, the truncated Green function can be Fourier 
transformed with respect to "fast' time t — t' , gaining 
a v slow ' time dependence and a linear correction in the 
oscillator velocity (See Eqs. (|A8IA9IA10|i ). 



1. Abiabatic Approximation: calculation of damping and 
fluctuating term 

We have now the tools to calculate the adiabatic cor- 
rections to the force acting on the mechanical oscilla- 
tor. For the sake of clarity, we rewrite here the adiabatic 
expansion of the Fourier transformed molecular Green 
function derived in the Appendix A fEq |A8p 



G r {u,t)~G r Q (u,t) + G{{u,t). 



(10) 



The explicit expressions of Gq and G\ are given by 
Eqs. (|A9IA10l) . Actually, in order to introduce dynam- 
ical effects on the force Eq.©, we have to calculate the 
adiabatic corrections to the lesser- Green function that is 
directly related to the occupation 



(n)(t) = -ihG < (t,t) 



E 

a=L,R 



h 2 T r 



|±/ Q H|G>,;)| 2 . 



(ii) 



From Eq. (|A9[) . one obtains, at zero-order, an expres- 
sion for the occupation of the same form of the static 
limit (Eq.(J7])) with the substitution E g o E g (t), acquir- 
ing a weak time dependence through the slow variable 
t. Adding the first order correction Eq. (|A10p into the 
Eq. (fTTj) , and neglecting terms proportional to the square 
velocity of the oscillator, we obtain 

(n)(t)~ h2r * I + 

a=L,R J 77 \ 



(12) 



In the end, the force given by Eq.© modifies to 

F(x) i-> F'(x, v) = F(x) ~ A(x)v, (13) 

where v = x is the velocity of the oscillator. We can con- 
clude that, both in equilibrium and in out-of-equilibrium 
conditions, the interaction with the leads introduces a 
dissipative correction term to the oscillator dynamics. 

Now, we observe that the introduction of a dissipa- 
tive term cannot be the unique dynamical effect for a 
classical dynamical system in contact with an environ- 
ment. It is well known that a fluctuating term should 
be included to take correctly into account the effect of 



the bath degree's of freedom. In our case, in order to 
include completely the effect of the " fast ' electronic en- 
vironment on the oscillator motion, we propose to take 
into account the fluctuations of the forced acting on the 
oscillator. These are induced by the intrisic "quantum' 
fluctuations of the electronic subsystem. 

We add to the average force contribution Eq.©, suit- 
able corrected by the damping term, Eq. flU)) . a stochas- 
tic fluctuating term able to take into account the effect 
of the electronic quantum fluctuations on the classical 
dynamics of the oscillator. Indeed, because we are con- 
sidering zero-temperature leads, one expects that these 
fluctuations are triggered by a finite bias voltage applied 
to the junction. As we shall see in the next sections, our 
approach is very accurate for junctions driven by a finite 
bias voltages, but it is far less accurate to describe the 
physics in the small bias (quantum) regime. 

We estimate the noise strength evaluating the average 
of the square fluctuation of the force over the electronic 
steady state. This fluctuating term is directly related to 
the fluctuation of the electron occupation 



(SF(t)8F(t')) = \ 2 {8h{t)8h{t')) 
\ 2 ({h{t)h{t'))-{h{t)){h{t')) 



(14) 



Decoupling the term (n(i)n(t')) with the Wick theorem, 
one obtains 

(SF(t)SF(t')) = X 2 h 2 G<(t' - t)G>(t - t') = 
= X 2 h 2 G<(t' -t)G>(t-t'), (15) 

where we have used zero-order time-dependent Green 
functions (as in Eq. (|A9I0 in order to take only first order 
corrections in the adiabatic ratio ^ . At this level of ap- 
proximation, we have obtained a multiplicative coloured 
noised According to the adiabatic approximation, we 
can further simplify the fluctuating term retaining only 
the zero-frequency component of the noise 



lim \ 2 K 



2 ' cfee i£(t -* 
D(x)8(t-t'), 



^G<(a, + £ )G>M 



(16) 



corresponding to electronic times scales comparable with 
that of the oscillator. In this way, one obtains a multi- 
plicative white noise term in the equation of motion for 
the oscillator. — The resulting Langevin equation for the 
oscillator dynamics becomes 



mx + A(x)x = F{x) + y/D(x)£(t), (17) 
(€(*))= 0, <£(*)£(*')}= 

where is a standard white noise term. Explicitly, the 
damping term A(x) is given by (from Eq. (|12p 'l 



A(x) 



AmujQ hug E p 

7T ~KThr 



E 

a=L,-Ft. 



2 + l] 2 



(18) 
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while the fluctuating term is (from Eq. (fT6 



D(x) 



muioEp hujQ 
n ~KY 



=L,-R \ 



arctan( 



E„ — Xx . 



KY 



-E„-\x 



nr 



,-Eg-Xa 
HT 



1], 



(19) 



where it is understood that ^2 a=L _ R K(a) — K(L) — 
K(R), for a generic function K(a). We note that 
Eqs. (fTT|) . (TTS)) and (ITO1) are identical to that obtained 
in Ref.23 — 25. Introducing a natural temporal unit 
to = 1/luq, the dimensionless damping A(x) and fluctuat- 
ing D(x) coefficients result proportional to the adiabatic 
ratio wo/r. As concerns the spatial dependence of the 
damping term, one can note in FigfT] (panels C-D) that 
it is almost localized on the position of the local min- 
ima of the static potential. The fluctuating coefficient, 
as shown in FigfTJ (solid (black) line in panels E-F), van- 
ishes at equilibrium (bias voltage Vuas = 0). Only ap- 
plying finite bias it becomes different from zero. In Figfl] 
(dashed (red) and dashed dotted (green) lines in panels 
E-F), we show that its spatial extension increases as the 
bias increases. Finally, one can note that A(x) and D(x) 
are almost independent by the ratio E g /E p . 



C. Numerical solution of Langevin equation: 
electronic observables and limits of the stochastic 
approach 

From the Langevin equation Eq. (fTT]) it is possible to 
derive the distribution probabilities P{x) and P(v) for 
the position and velocities variables of the oscillator. We 
have evaluated them solving the second order stochastic 
differential equation with the 4-order stochastic Runge- 
Kutta algorithm developed by R. L. Honeycutti 35 i 36 First 
of all, as suggested in Ref.37, in order to solve our second 
order ordinary differential equation with multiplicative 
white noise, we decompose the problem in a set of three 
first order differential equations. The third one takes 
into account the effect of spatial dependence of the noise 
involving a non multiplicative noise term. For our sim- 
ulations we have fixed a time step t s = O.lr (r = 1/wo) 
and set long simulation times up to T = lO 9 ^. Within 
these settings, the algorithm shows an excellent stabil- 
ity in the whole range of model parameters. In order to 
construct our histograms, we have sampled the values of 
x(t) and v(t) every 100 time steps. We have therefore 
obtained the distribution probabilities for the stationary 
state of the oscillator. 

Given our assumption about the separation between 
the slow ionic (vibrational) and fast electronic (tunnel- 
ing) timescales, the problem of evaluating a generic ob- 
servable (electronic or not) of the system reduces to the 
evaluation of that quantity for a fixed position x and ve- 
locity v of the oscillator, with the consequent averaging 
over the stationary probability distributions, P(x) and 



P(v). Therefore, for a generic observable which depends 
only on position, 0(x), the averaged quantity is 



(0(x)) = / dxP(x)0{x), 



(20) 



while, for an observable which depends on velocity vari- 
able only, one has 



(0(v)} = / dvP(v)0(v). 



(21) 



In our case, the current, the spectral function and the 
electronic occupation depend only on the position vari- 
able 



I(x) 
A Spec {uj,x) 



nr 



A S pec(uJ,x), (22) 



(hcu- E g - Xx) 2 + h 2 Y 2 /^ 

Ha — E g — XX 



(23) 



(n)(x) = - + — arctan 



2?r 



a=R,L 



KY/2 



(24) 



The position distribution probabilities P{x) have been 
already discussed by authors of Ref.24 in the extremely 
strong coupling regime E p » hY » huj a . They ana- 
lyze the case where the static potential shows two sym- 
metric or asymmetric wells separated by a very high bar- 
rier. In this regime, solving numerically the Fokker-Plank 
equation of the problem, they estimate the switching- 
rates by evaluating the escape times from each well of 
the generalized potential. Indeed, this point is inter- 
esting for clarifying the role of electron-phonon interac- 
tion in the appearance of a bistable behaviour in single 
molecule tunneling devices. One of the results of this 
paper is that the multistability and hysteretic behaviour 
in the current- voltage characteristic disappear if the dy- 
namical effects of the oscillator motion are taken into 
account. To clarify this point, we focus here on the case 
(already considered in Ref.24) where the switching times 
between different oscillator potential wells are very long, 
and the oscillator jumps between two states (see panel A 
of FigfS]) corresponding to very small electronic currents. 
In order to explore the same regime of parameters, in 
our approach very long simulation times as T = 10 9 t s 
are necessary for sampling the entire phase space expe- 
rienced during the dynamics. Nevertheless, as shown in 
Panel B of Fig[2J we get an excellent agreement with 
Pistolesi et al. results. It is interesting to note that in 
the paper of Ref.24, the authors consider a broadening Y 
which is twice our values (we show in the caption of FigfS] 
the comparison between the simulations taking correctly 
into account this factor). In the small bias regime, we ob- 
serve a strong suppression of the current. The oscillator 
spends a long time in each potential well, suddenly jumps 
into the other and then comes back in the same way (see 
panel A in Figf2]) . For clarity, we show in panel C of Figf2] 
the corresponding position distribution probability P{x). 
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FIG. 2. (Color online). Panel A: Solution of the Langevin 
equation Eq. (|17p in the extremely strong coupling regime 
E p » hT » huj (hw /2E p = 10~ 3 ) for hr/2E p = 0.08, 
E g = Ep and eVBias/2E p = 0.1. Panel B: Current (er units) 
voltage (eVbias in 2E P units) characteristic for the same value 
of hV/2E p as above. Solid (red) curve is drown from Ref.24, 
square line indicates our dynamic simulation and dashed line 
indicates static I-V. Panel C: Dimensionless position distribu- 
tion probability for the same values of parameters as in Panel 
A. The dimensionless position x, time i and distribution func- 
tion P are defined as x = x/xo, i = t/to, P = P/(l/xo), with 
xq — and t = l/wo, respectively. 



FIG. 3. (Color online). Panel A: log-plot of dimensionless 
velocity probability distribution function vs. v 2 , at different 
adiabatic ratios (the values of ujq shown in the figure are in 
P units), fixed bias voltage Vwas = 0.1 and different gate 
voltages and EOC strengths (not shown in the graph). The 
dotted (red) lines indicate that curves have a good linear fit. 
Panels B — C — D: log-plot of velocity probability distribution 
function vs.v 2 for Vbias = 0.1, Vbias = 1-1, Vbias = 2.1, respec- 
tively. The dashed (red) line indicates linear fitting. Dotted 
(green) and dash dotted (blue) lines indicate polynomial fit- 
ting of 2nd and 4th degree. Vuas values are expressed in hT/e 
units. The dimensionless distribution function is defined as 
P = P/(mujo/X), while v 2 is expressed in (X/mujo) 2 units. 



The maxima of P(x) correspond to two small current car- 
rying states: the position of the molecular energy level is 
far above (E g ~ E p ) or below (E g ~ — E p ) the chemical 
potential of the leads. For sufficiently large bias voltage, 
as discussed in subsection A, appears a third minimum 
in the static potential. This minimum corresponds to 
a large-current carrying state determining a continuous 
enhancement of the current, against the abrupt discon- 
tinuity (or hysteresis) which would been obtained in the 
static approximation (dashed line in panel B of FigfJJ. 



1. Non gaussian features of P(v) and study of the average 
kinetic energy of the oscillator 

In this section we focus our attention on the oscillator 
observables 0(v) which depend on the velocity v. We re- 
mark that the oscillator is coupled to the electronic bath 
only through the interaction term Hi n t, Eq.([5]). As the 
bias voltage increases, this bath is strongly driven out 
of equilibrium. It is therefore important to analyze the 
effect of the electronic subsystem on the oscillator dis- 
tribution probability P(v) as a function of the bias volt- 
age. In the small bias regime, regardless the value of the 



gate voltages E g and the coupling E p , as shown in FigJ3] 
(panel A) for different adiabatic ratios (from ujo = 10~ 3 
to wo = 0.25), the velocity distribution probabilities P(v) 
are gaussian. In this regime, the nonequilibrium elec- 
tronic bath behaves like a conventional bath for the os- 
cillator with an " effective ' temperature linearly propor- 
tional to the bias voltage. As described in the inset of 
FigHl at arbitrary E p and gate voltages the kinetic en- 
ergy curves show a common linear trend at small bias 
with a slope V5,j as /4 in agreement with Mozyrsky et al. 
(we get Vuas/ 8 because we choose a broadening HT half 
that used in Ref.23). As we increase the bias voltage, 
the (logP(u) vs. v 2 ) plot starts to deviate from a linear 
trend, as shown in FigOl panels B — C — D. This be- 
haviour indicates that the oscillator dynamics cannot be 
simply reduced to an effective temperature in this regime, 
pointing to a very significant role of the dynamical effects. 

In the adiabatic approximation, the average kinetic 
energy of the oscillator has an important role. It de- 
scribes the effect of the "back-action' of the nonequi- 
librium electronic bath on the oscillator dynamics and 
can be used, as shown below, as a tool to assess the 
validity of the adiabatic approximation. We show in 
FigH]the behaviour of the kinetic energy {Ejci n } for dif- 
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FIG. 4. (Color online). Main: plot of average kinetic energy 
{Exin) as function of the bias voltage at fixed adiabatic ratio 
uio/T = 0.1 and gate voltage E g — 0, for different interaction 
strengths E p : E v — 0.1 square (black) curve, E p — 1.0 circle 
(red) curve, E v = 2.0 triangle (green) curve, E v = 3.0 star 
(blue) curve. Two constant energy lines E = Tiujq/2TiT — 
0.05 (dashed) and E = hT — 1 (dotted) are also plotted. 
Inset: Average kinetic energy (Eicin) for small bias voltages 
for the same parameter values of the main plot. The dotted 
(magenta) line indicates the linear approximation eVbi as /8 
derived in Ref.23 (we choose a broadening hT half that used 
in Ref.23). All the quantities ({Exin}, E p , E g and eVw as ) are 
in unit hT. 

ferent interaction strengths E p as function of the bias 
voltage. First of all, we note that, regardless the val- 
ues of E p , for Vbias — all kinetic energy curves show 
(Ekin) — 0. At equilibrium, we can say that the oscillator 
v thcrmalizes ' to the temperature of the electronic bath 
(T e i = 0). We have also plotted two constant energy lines 
that specify the range of validity of our approximation, 
E = hojo/2 <~ and E = fiT. At intermediate bias 

values, the curves show a departure from the common 
linear behaviour observed in the small bias regime, more 
evident as the interaction strength increases. The kinetic 
energy curves corresponding to E p — 2.0 and E p = 3.0 
show an interesting plateau at intermediate bias where 
increasing the bias does not produce an increase of the 
average kinetic energy. Actually, at E p = 3.0, we find 
even a very slight decrease. We also note, in the same 
regime, that the velocity distribution probabilities are 
not gaussian. 



2. Limits of the adiabatic approach 

As mentioned, we can use the average kinetic en- 
ergy of the oscillator to fix the range of validity of the 



FIG. 5. (Color online). Diagram for the range of validity of 
classical approximation at fixed adiabatic ratio uiq/T — 0.05 
(the value of ojq shown in the figure is in F units). The dashed 
(black) line indicates the QR-CAR crossover for E g = and 
E g = 1. The dotted (red) and dashed dotted (green) lines 
indicate the CAR-CNAR crossover for E g = and E g = 1, 
respectively. E p , E g and eVtias are expressed in unit KT. 



adiabatic approximation. If this energy is lower than 
the characteristic Debye temperature of the oscillator 
(Exin) < Tiojq/2 ~ we actually explore a region, 

as discussed in Ref.40, where quantum correlation effects 
can not be disregarded. We call this region Non Classical 
or Quantum Region (QR). If k B T D /2 < (E Kin ) < HT, 
that is the kinetic energy is lesser than the characteristic 
energy scale of the electronic degrees of freedom and si- 
multaneously greater than characteristic Debye tempera- 
ture, a huge number of vibrational quanta (phonons) are 
excited in the system. We call this region Classical Adi- 
abatic (CAR). When the average kinetic energy of the 
oscillator exceeds the characteristic energy scale of the 
electron dynamics (EKin) > frT, we clearly are going be- 
yond the limit of adiabatic approximation we start with. 
We define this region Classical Non Adiabatic (CNAR). 
We expect that in the CAR our approximation is very 
accurate. By using this data, we are now able to build 
up a diagram for the validity of classical approximation 
in the plane (E p -Vbias) for different values of gate volt- 
ages (Fig[5]) and different adiabatic ratios (Figj6]) . It is 
interesting to note that, in FigJSJ the QR-CAR crossover 
line is almost independent from the gate voltage in the 
limit of small adiabatic ratio. Instead, the CAR-CNAR 
crossover line is slightly dependent from the gate voltage 
showing an enlargement of the CAR with E g . Globally 
we note that, apart for the QR (small bias), the CAR 
enlarges as one increases the electron oscillator coupling. 

As expected, as we increase the adiabatic ratio, the 
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FIG. 6. (Color online). Diagram for the validity of clas- 
sical approximation at fixed gate voltage E g = (asym- 
metric static potential) for different adiabatic ratios ujo/T = 
0.01 — 0.05 — 0.1 — 0.25 (the values of ujq shown in the figure 
are in T units). E p , E g and eVuas are in TiT units. 



QR expands reaching great values of bias voltage 14j as , 
FiglHl In particular, observing the down-triangle (blue) 
curve (adiabatic ratio uj^/T = 0.25), one notes that the 
QR-CAR crossover line reaches a * v maximum ' ' in cor- 
respondence of E p ~ 2 and Vuas — 2.8. For E p < 2, the 
bias values identifying the QR-CAR crossover increase 
as the voltage increases. For couplings E p > 2, we note 
an inversion of this behaviour: the CAR starts to extend 
for a very large area of the diagram except for a narrow 
region at small bias (QR) and for a region at bigger bias 
values (CNAR). This means that, even for intermediate 
adiabatic ratios, we need sufficiently strong couplings E p 
in order to obtain a predominant CAR in the validity 
diagram. Moreover, this is due to the fact that the node 
between kinetic energy curves and the Debye line occurs 
in the non monotonic intermediate bias region (see FigSJ) . 
On the other hand, the CAR-CNAR crossover line is al- 
most independent from the adiabatic ratio (for not too 
large interaction strength). This is what we expect from 
physical grounds and constitutes a self-consistent check 
of our approximation. 



3. Electronic transport properties 

We can now analyze the electronic transport prop- 
erties resulting from the average over the dynamical 
fluctuations of the oscillator motion. We first study 
the conductance-voltage curves as function of the EOC 
strength (FiglTJ), then we show how the dynamical fluc- 
tuations strongly renormalize the infinite mass approxi- 




bias 



FIG. 7. (Color online). Panel A: Conductance (in units Go = 
fr) in the static approximation as function of bias voltages, 
for ojo/F = 0.05, Eg = and different interaction strengths 
E p = 0.05,0.5,1.0,2.0,3.0. Panel B: Dynamical correction 
to the conductance for the same parameter values of panel 
A. The value of ujq shown in the figure is in T units while 
all other quantities (E g , E p and eVbias) are expressed in TiT 
units. 



mation results studying the I-V curves for different adi- 
abatic ratios (Fig|5J|. Finally, we investigate the depen- 
dence of the kinetic energy and of the I-V characteristic 
as function of gate voltage studying the properties of the 
junction as function of an electric field orthogonal to the 
source-drain direction (FigJSJ). 

In Fig[7]we show several conductance curves for differ- 
ent interaction strengths, E p = 0.05 — 3, at luo/T = 0.05 
and Eg = 0. The comparison between static (panel A) 
and dynamical (panel B) approximation is very interest- 
ing. The static solution shifts the non interacting reso- 
nance by a quantity proportional to the polaronic energy 
E p (panel A). As one can see, this effect strongly reduces 
the small bias conductance. The dynamical correction, 
on the other hand, reduces the polaronic shift compared 
to the static curves and also broadens (as a result of 
the very broad nonequilibrium distribution probabilities 
P{x)) the electronic resonance. In the intermediate bias 
regime, we note a strong enhancement of the conduction 
far from the electronic resonance where a very small cur- 
rent is observed in the static approximation. Moreover, 
including the dynamical fluctuations, the reduction of the 
small bias conductance is less pronounced. We note also 
that our dynamical approximation is close to the static 
solution in the small bias regime, while is substantially 
different in intermediate one. The dynamical corrections 
strongly renormalize the static results even for small adi- 
abatic ratios. 
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FIG. 8. (Color online). Panel A: Electronic occupation as 
function of bias voltages, for E g = 0, E p = 2 and different 
adiabatic ratios uio/T — 0.05,0.25,0.5, 1.0. Panel B: Current 
voltage characteristic for the same value of the parameters of 
panel A. The values of wo shown in the figure are in V units 
while all other quantities (E g , E p and eVbias) are expressed 
in hT units. 



FIG. 9. (Color online). Panel A: Current as function of gate 
voltages, for loq/T — 0.05, E v = 0.25 and different bias volt- 
ages Vbias = 0.1,1.0,2.0. Panel B: Plot of average kinetic 
energy as function of the gate for Vbias = 0.1,1.0,2.0. The 
value of ujq shown in the figure is in T units, while E g and E p 
are expressed in HT units. Vbias is expressed in hT/e units. 



We analyze in FigJ5] the behaviour of the electronic 
occupation (panel A) and current voltage characteristic 
(panel B) at strong coupling E p — 2, for different adia- 
batic ratios uj /T = 0.01, 0.1, 0.25 and at E g = 0. In the 
small bias regime, as a result of strong electron-oscillator 
interaction, the molecular level renormalizes itself far be- 
low the chemical potential of the leads. We note a large 
difference between the non interacting occupation value 
((h) ~ 0.5) and the interacting one ((n) ~ 1). As one 
increases the bias voltage, many charges are pumped out 
the molecular v dot ' . In the large bias regime the sta- 
tionary charge quantity in the molecular " dot ' reduces 
approaching the non interacting value ((h) ~ 0.5). The 
nonequilibrium broadening of the distribution probabili- 
ties P(x), then, induces a strong reduction of the conduc- 
tion threshold with respect to the static solution (solid 
magenta curve in Fig|8] panel B). We note a small varia- 
tion of physical properties with respect to the adiabatic 
ratio at intermediate voltages, in the CAR in correspon- 
dence to non-gaussian regime of the distribution proba- 
bilities. 

In many molecular transport experiments, one records 
the current or the conductance varying an electric field 
applied on the molecule (orthogonal to the source-drain 
direction) at fixed source-drain voltage. In the panel 
A of Fig[9l it is shown the current for different bias 
voltages at moderately small electron-oscillator coupling 
(E p = 0.25) as function of gate voltage. In this regime 
we have no bistability in the model (E p < KT). We note 



that the static and dynamical approximations agree in 
the small bias regime (solid and square lines in panel A 
of Fig[5]). Increasing the bias voltage, the dynamical cor- 
rection becomes more important showing a suppression 
of the current for small E g . This effect is caused by the 
spectral weight broadening due to the average over the 
position distribution probabilities P(x). From panel B 
of FiglHl (square (red) line), we learn that, in the small 
bias regime, the kinetic energy is independent of the gate 
voltage, while, as the voltage increases, it shows a sym- 
metric drop with respect to polaronic energy E g = E p , 
corresponding to the symmetric regime. The I-V curve 
also shares this symmetry. This effect can be explained 
observing that, when the s bare' molecular level and the 
rcnormalized one are both in the bias window, the energy 
associated to the electronic current flow is more efficiently 
exchanged with the oscillator. When the electronic res- 
onance is far above or below the chemical potential of 
the leads there is a less effective coupling between the 
oscillator and the electronic subsystem. 



III. THE TWO-SITE SSH MODEL 

The first step towards a more realistic description of 
a molecular junction is to consider a model Hamiltonian 
composed by two sites connected by an internal hopping 
t. In particular dimer molecules,^ this hopping can be 
controlled by a vibrational mode which assists the elec- 
tron tunneling through the two molecular sites. In this 
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case, a guess for the molecular Hamiltonian is given by 
tffSff = Eg {S[di + d\d 2 ) - t{x)(d\d 2 + h.c), 

(25) 

where we consider, as in the SSH model, an electron hop- 
ping 



t(x) = t — ax 



(26) 



depending linearly on the lattice displacement x asso- 
ciated with the intermolecular vibrational mode. The 
molecular sites have a common energy E g and are de- 
scribed in terms of creation (annihilation) operators 



d}(di 



1,2. The SSH model was indroduced to 



describe the transport properties of conducting polymers 
(e.g. polyacetylene22) and the two site case represents 
the shortest version of a molecular wire^ A generaliza- 
tion of this two site model was proposed in Ref.32 for the 
study the electron transport of dimer molecules interact- 
ing with a single internal vibrational mode. 

Most molecular devices studied experimentally so 
fa r 2 ' 3 ' 38 have been weakly coupled to the leads. This cor- 
responds to the bare tunnel broadening hT of molecular 
electronic levels smaller that the energy required to ex- 
cite one oscillator quantum (phonon) Hujq . In the strong- 
coupling regime, when the electron-oscillator interaction 
energy E p exceeds Hujq, the physics is governed by the 
Franck-Condon effect>i2~— , i.e. the tunneling of an elec- 
tron onto the molecule with the simultaneous emission 
or absorption of several phonons is more probable than 
elastic tunneling. The current as the function of voltage 
exhibits steps separated by huia/e^ and the conductance 
shows phonon sidebands^ 

As in the AH model, we study here the case of slow 
phonons, loq << T, coupled to a molecular junction 
driven by a finite bias, in particular for eVuas > huj^. 
As further approximation, we consider the dynamics of 
the vibrational mode " " classical ' ' . 

The structure of the SSH model is very interesting. 
The direct coupling of the electron-oscillator interaction 
to the intermolecular hopping t suggests that the role of 
the dynamical fluctuations becomes crucial to determine 
the physical scenario. The total Hamiltonian is 



HtOT = Hel-SSH + H D 



where 



el-SSH 



H 



SSH 
Mol 



Tun 



(27) 



(28) 



with Hieads and H osc given by Eq.Q and Eq.Q, respec- 
tively. The tunneling Hamiltonian H-Tun is given by 

H Tun = ^(Vk,LC ktL di + h.c.) +^2(Vk,Rcl R d2 + h.c), 



k,L 



k,R 



(29) 

indicating that the left (right) lead is coupled only to the 
molecular site 1(2). In real space, the molecular Hamil- 
tonian Hfjjj 1 is not diagonal. We therefore perform a 
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FIG. 10. (Color online). Sketch of junction within the SSH 
model in an energy scale. 



transformation which diagonalizes the molecular isolated 
problem 



Si 



72 



d\ + 4 
~7T~ 



V2 



(30) 



with the same transformation for corresponding anni- 
hilation operators. This transformation leaves invariant 
Hieads but changes Hffjj 1 and Hrun- Explicitly we have 



HmoI ~ e 7i i x ) c \i c 7i ^~ £ 72 ( x )c\ 2 ^72 ' 



(31) 



H 



Tun 



E 



i=l,2 



V k , L 4 



V2 



Vk,R 4 -ii, \ 

^2 C k,R C li ' n - C -> 



where 



(x) = e + (t — ax), 



,0*0 



(t — ax). 



(32) 



(33) 



As one can see, the above transformation allows us to 
take into account exactly the intermolecular hopping's 
effect. The molecular Hamiltonian H^f Q f (Eq.(l5T1l) is 
equivalent, at fixed x, to that of a non interacting two 
level system. In FigfTUl a schematic picture of the junc- 
tion in an energy representation is shown. We observe 
that there are two electronic resonances, corresponding 
to a v bonding ' and " anti-bonding ' states whose position 
is renormalized by the electron-oscillator interaction. 

From now on we work in the energy space for conve- 
nience. 
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In the following, we will (subsection A) first ana- 
lyze the coupled electron-oscillator problem within the 
SSH model in the limit of infinite mass for the oscilla- 
tor. Then, we will construct, as done in AH model, the 
stochastic Langevin equation for the dynamics of the os- 
cillator. In the subsection B we will describe the numer- 
ical results. 



A. Out of equilibrium Born-Oppenheimer 
approximation: infinite mass (static) case 

As in AH model, performing the limit m t— > oo, at 
zero-order static approximation, we neglect the kinetic 
energy of the oscillator. The electronic dynamics, with 
the oscillator displacement free parameter, is there- 

fore equivalent in the energy space to a non-interacting 
two level problem with energy levels renormalized by the 
"polaronic' shift —ax, Eq. (|33[) . In what follows, we con- 
sider the case of symmetric coupling of the molecule to 
the leads hT^ = HYr in the wide-band approximation. 
Here, we briefly show how to calculate the generalized 
potential of the oscillator coupled to the double " dot ' 
molecular junction. 

Within the Keldysh formalism, we use the equation of 
motion approach to calculate the molecular Green func- 
tions in stationary nonequilibrium conditions. In the zero 
order static approximation, we have the following equa- 
tion of motion for the molecular retarded Green function 



ih-Sr - e 11 (x) + i 

1 A 



4 



, nr L -nr r 



^2,l(*>*) Gvj, 2 (M') 



= *(*-0(j ?), (34) 



which acquires a 2 x 2 matrix structure. A similar equa- 
tion is valid also for the advanced Green function. In 
the hypothesis of symmetric coupling with the leads, one 
obtains two separate problems for the molecular energy 
levels e 11 (x) and s l2 (x), respectively. The diagonal el- 
ements of the retarded Green function in Fourier space 
are 



1 



fiui — e~ H (x) + ^^ 



hr L +nr n 

4 



i = 1,2 



(35) 



while the non diagonal terms are zero. 

The lesser matrix Green function is instead given by 

G < {w,x) = 



4 ^ K ~ "«) G 2,2 G ?,i (n L +n R )\G r 2j2 \ i ) ^> 

where, for sake of simplicity, we have dropped the fre- 
quency uj and the displacement x dependence. The di- 
agonal terms of the lesser Green function are directly 



related to the electron " " densities ' ' (these obviously not 
correspond to the densities in real space) 



)(x) = - + — arctan 
/w 2 2vr ^ 



a=R,L 



Mo 



£ 7i (x) 



TiT/A 



1,2. 



(37) 



For sake of clarity, we show here that the population in 
real space of the left and right molecular sites are ex- 
pressed in terms of lesser Green functions fEq. flMl) ') 

(ni){x) = \j ^(G^ + G^ + i-iy^iG^ + G^] 

(38) 

where i = 1,2 (site 1 is the left site, site 2 the right one). 
The force exerted on the oscillator is given by 



Fssh = ~kx + a((n. 



71/ 



'72 



»(*)■ 



(39) 



Taking care of Eq. (|37j) and Eq.([39|). one can straight- 
forwardly compute the expression of the generalized po- 
tential in nonequilibrium conditions (hr = — eVbi as /2, 
fi L = eVbias/2) 



Vssh(x) = -kx 2 



arctan 



'Mo 
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(40) 



This generalized oscillator potential depends parametri- 
cally on the new electronic energy scale introduced in 
the problem: the intermolecular hopping t " " hidden 
in e 7i (x), see Eq. (|3"3"l) . Furthermore, it depends on the 
polaron energy, E p , the gate voltage E g , and the bias 

In Fig llll we present some features of the generalized 
potential Vssh(x) which will allow us to understand the 
effect of the nonequilibrium electronic system on the 

static ' ' stretching of the oscillator (solutions of the 
equation Fssh = 0). Moreover, this will help us to clar- 
ify the role of the dynamical effects in the transport prop- 
erties that we will show later. 

We focus here on the weak coupling (E p /frT « 1) 
regime where moreover the intermolecular hopping t is 
larger than the coupling TiT of the molecule with the 
leads. In the panel A we show the generalized poten- 
tial of the SSH model at fixed EOC strength, E p = 0.2, 
intermolecular hopping t = 2.0, as function of the bias 
voltage Vbias- One can observe that, as the bias increases, 
the position of generalized potential minimum goes from 
x ~ — 1 (corresponding to (n 71 ) — (n 72 ) ~ — 1) to x ~ 
(corresponding to (n 7l ) — (n 72 ) ~ 0). The oscillator 
switches from a full stretching configuration (x ~ — 1) 
to a no-stretching one (x ~ 0). At equilibrium (solid 
(black) curve of Panel A of Fig lTTj) . we have a physical 
situation where the renormalized anti-bonding electron 
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FIG. 11. (Color online). Panel A: Spatial dependence of the 
dimensionless generalized static potential Vssh(x) at u)q/T — 
0.1, E p = 0.2, gate voltage E g = 0, intermolecular hopping 
t = 2.0, for different values of the bias voltage: Vbias = 0.0 
(solid (black) curve), Vbias = 3.5 (dashed (red) curve), Vbias = 
4.0 (dotted (green) curve), Vbias = 6.0 (dashed dotted (blue) 
curve). The vertical lines indicate the position of the minima 
of the potential. Panel B: Same as above for cjo/F = 0.1, E p — 
1.4, t = 0.2, gate voltage E g — 2 and different values of the 
bias voltage: Vu as = (solid (black) curve), Vuas = 4 (dashed 
(red) curve), Vbias = 8 (dotted (green) curve). The potential 
is expressed in TiT units (Vssh = Vssh/KT). Vbias values 
are expressed in HT/e units, where e is the electron charge. 
Eg, E p and t are expressed in TiY units. The dimensionless 
position variable x is defined as x — x/xq with xq = — 



level e lx (x) is above the chemical potential of both leads, 
while the bonding one s 12 (x) is below them. The clas- 
sical ""spring'' is fully compressed (x ~ —1) and this 
corresponds in real space to molecular sites half-filled 
((ni) ~ (712) — 0.5). Studying the electronic popula- 
tions of left (1) and right (2) molecular sites (Eq.(|38jl). 
one can observe that, if we increase the bias voltage, the 
left site starts to empty, while the right one populates, 
reaching, for eV b * as /2 ~ t ~ ax(V b * as ) (hopping value 
properly renormalized) , a small difference of population 
roughly equal to (ni) — (n 2 ) — —0.1. For sufficiently 
large bias, the molecular level populations tend again to 
the common value 0.5. As we shall see in next section, 
the inclusion of the dynamical effects allows to clarify 
the physical picture arising from the above description, 
in terms of an energy balance between the electronic and 
oscillator subsystems. 

At static level, it is also interesting to discuss the 
extremely strong coupling regime E p > hT > t, for 
gate voltage E g = 2.0 (panel B). In this case, at equi- 
librium, we are describing a physical situation where 



the renormalized bonding and anti-bonding electron lev- 
els are both above the chemical potential (Vbias = 0). 
The molecular sites in real space are both almost empty 
((ni) ~ (712) — 0), and the oscillator is in a no-stretching 
configuration x ~ (solid (black) curve). Increasing 
the bias voltage, the generalized potential develops dif- 
ferent minima. At intermediate bias, one can observe 
two asymmetric minima near x ~ —0.5 and x ~ 0.5, 
separated by a potential barrier. In this regime, the 
minimum corresponding to x ~ —0.5 prevails (dashed 
(red) curve) and the non-interacting real space popula- 
tions (111) and (n2) are asymmetrically distributed be- 
tween the two sites ((n\) ~ 0.8, (n 2 ) ~ 0.2). In- 
stead, the interacting real space populations have the 
same value (m) ~ (n 2 ) — 0.25, corresponding to a very 
large current-carrying configuration. In the large bias 
regime only the minimum x ~ corresponding to a 
small-current-carrying configuration survives. Including 
the interaction effects, the left site results almost filled 
(rii) — 0.9, while the right one almost empty (n 2 ) — 0.1, 
showing that, as result of the strong electron-phonon in- 
teraction, the bias voltage does not manage to deplete 
both molecular sites. As we shall see later, the features 
of the static potential obtained in this case determine the 
possibility to observe in the I-V a strong Negative Dif- 
ferential Resistance, when the dynamical effects of the 
oscillator are neglected. 



B. Adiabatic Approximation: calculation of 
damping and fluctuating term 

As we have discussed after the Eq. (IMl) , the assumption 
of symmetric coupling to the leads allows to disentangle 
in the energy space the problem for the molecular bond- 
ing and anti-bonding levels e 7i (x) . Repeating site- by-site 
the construction introduced in the previous section for 
AH model, we can straightforwardly set for our two site 
SSH model a Langevin equation for the oscillator dynam- 
ics, very similar to that derived in AH model. The new 
coefficients, F(x), A(x) and D(x) are given by 



F(x) = -kx + X±- J2 E 

a=R,Li=l,2 
a=L,-Ri=l,2 
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where in the first line of Eq. (|43l) we have dropped the 
frequency u and the displacement x dependence in the 
Green functions. We end this section briefly discussing 
some of the peculiarities of the damping function A(x) 
and of the fluctuating term D(x). As regards the damp- 
ing term (panel A and C in Fig. (fT2|) ). one can observe 
that is located in suitable points and is strongly space 
dependent. It is interesting to note that, as in the AH 
model case, it survives also for Vuas = (solid (black) 
curves in Panels A-C). In this case, one can also note 
that, for E p << t (panel A), A(x) is almost zero in the 
interval mostly explored in the dynamics (—2 < x < 2), 
while, for E p » t (panel C), shows two pronounced 
peaks in that interval. As one can see, the position of 
A(xYs maxima is strongly bias dependent. 

As concerns the fluctuating term (Panels B and D in 
Fig. (TT2"10 . one can note that, as in the AH model, it is 
identically zero at equilibrium (Vhias = 0). When the bias 
increases, it becomes almost different from zero in the 
region mostly explored in the dynamics (—2 < x < 2). 
In Panel B of Fig[l2] one can observe that the spatial 
extension of the fluctuating term increases as the bias in- 
creases, while, in Panel D, in the interesting strong cou- 
pling regime (E p » t), it shows a maximum for x = 0, 
the no-stretching equilibrium state of the oscillator. It is 
important to stress again here that the space dependence 
of these terms determines the non-gaussian character of 
the distribution probabilities P(x) and P(v) of the oscil- 
lator. 




FIG. 12. (Color online). Panels A-B: Spatial dependence 
of the dimensionless friction coefficient A(x) and fluctuating 
term D(x) at ljo/F = 0.1, E p = 0.2, gate voltage E g — 0, 
intermolecular hopping t = 2.0, for different values of the bias 
voltage: Vbias = 0.0 (solid (black) curve), Vbias = 3.5 (dashed 
(red) curve), Vbias = 4.0 (dotted (green) curve), Vbias = 6.0 
(dashed dotted (blue) curve). Panels C-D: Same as above for 



wo/r = o.i, 



1.4, t — 0.2, gate voltage E g — 2 and 



different values of the bias voltage: Vbias = (solid (black) 
curve), Vbias = 4 (dashed (red) curve), Vbias = 8 (dotted 
(green) curve). The friction coefficient is expressed in mcoo 
units (A — A/mu)o) while the fluctuating term in A 2 /ujq units, 
(D = D/(\ 2 /wo)). Vbias values are expressed in hV/e units, 
where e is the electron charge. E g , E p and t are expressed in 
TiT units. The dimensionless position variable x is defined as 
x — x/xn with Xn = — ^-r. 



C. Analysis of Numerical results 

As done for the AH model, we here show the results 
arising from the numerical simulation of the Langevin 
equation of the SSH model. We evaluate the fundamen- 
tal ingredients of the adiabatic approximation: the distri- 
bution probabilities for the oscillator. These allow us to 
calculate the dynamical properties of the oscillator (aver- 
age kinetic and potential energy) as well as the electronic 
transport properties of the molecular junction. 



1. Study of the average kinetic energy of the oscillator and 
limits of the Adiabatic approach 

First of all, we study the behaviour of velocity dis- 
tribution probabilities P(v) resulting from the solution 
of the Langevin equation associated to the SSH model. 
As in AH model, we have verified that in the small bias 
regime, regardless the value of the gate voltages E g , the 



electron-oscillator coupling E p and the hopping t, the di- 
mensionless velocity distribution probabilities P(v) are 
gaussian. The introduction of a new energy scale in the 
problem does not much modify the physical picture we 
obtained in AH model in the small bias regime: the non- 
equilibrium electronic bath behaves like a conventional 
bath for the oscillator with an effective temperature lin- 
early proportional to the bias voltage. In particular, in 
the SSH model case, it is worth noticing that the average 
kinetic energy exhibits a slope twice that found in the 
AH model. This is a consequence of the transformation 
Eq. (f3"0)) we have applied on the total Hamiltionian, that 
renormalized the tunneling amplitudes with the leads, 
Vk,a i-> Vk,a/v2- From the physical point of view, we 
find that the two electronic channels independently con- 
tribute to the oscillator effective temperature, showing 
that the problem is equivalent to the sum of two single- 
site junctions. 

As we increase the bias voltage, the (log(P(w)) vs. v 2 ) 
plot starts to deviate from a linear trend, so that, even 
in SSH case, the oscillator dynamics cannot be simply 
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FIG. 13. (Color online). Diagram for the validity of classical 
approximation at fixed adiabatic ratio wo/P = 0.1, E g = 0, 
for different values of intermolecular hopping t — 0.2 — 0.6 — 
1.0 — 2.0. The value of ujq shown in the figure is in F units, 
while all other quantities (E g , E p , t and eVuaa) are expressed 
in hV units. 



reduced to an effective temperature in the intermediate 
bias regime. 

We also find that, for < i/£ p < 1 and up to very 
large values of the bias voltage, the average kinetic en- 
ergy shows a behaviour qualitatively similar to that of 
AH model (FigBJ. In this regime, we can conclude that 
the dynamical fluctuations of the oscillator motion do not 
s see ' the double " dot ' structure of the electronic molec- 
ular junction. If t/E p >> 1, as we will discuss later, the 
average kinetic energy shows an interesting non mono- 
tonic behaviour in the intermediate bias regime (see be- 
low, Fie fM]) . 

The systematic analysis of the average kinetic energy 
allows us to build up a diagram for the validity of classi- 
cal approximation in the plane (E p -Vbias), as done for AH 
model (FigfTSI). In particular, we study the validity dia- 
gram for different values of intermolecular hopping t and 
at fixed adiabatic ratio and gate voltage. In this case, 
it is interesting to note that QR-CAR crossover line is 
almost independent by the intermolecular hopping in the 
limit of small adiabatic ratio. Joining together the results 
obtained for the AH validity diagrams (FigJS]and FiglH]), 
we can conclude that, in the limit of very small adiabatic 
ratios, the QR-CAR crossover line is completely indepen- 
dent by the other energy scales considered in the prob- 
lem. The CAR-CNAR crossover line is instead slightly 
dependent on t showing the expansion of the CAR. As the 
intermolecular hopping t increases, bigger values of bias 
voltage are needed to get average kinetic energy values 
greater than energy coupling to the leads, (Exm) > hT. 



FIG. 14. (Color online). Panel A: Average potential energy as 
function of the bias for different values of intermolecular hop- 
ping t — 0.2 — 0.6 — 1.0 — 2.0. Panel B: Average kinetic energy 
as function of the bias for values of intermolecular hopping 
as in panel A. We note that the introduction of new energy 
scale makes the overall energy (E) a decreasing function of 
the voltage for intermediate values. The value of loq shown 
in the figure is in Y units, while all other quantities (E g , E p , 
(Exin), {Epot}, ksT, t and eVuas) are expressed in TiY units. 



In this case, the intermolecular hopping t plays the same 
role as the gate in AH model (see, FigJS]). 

A new feature which was not observed in the AH model 
is the appearance of small QR for sufficiently small cou- 
pling E p , at intermediate bias voltages (Fig fTSl) . For 
strong enough electron-phonon coupling E p , these re- 
gions disappear. This feature can be understood ana- 
lyzing the behaviour the average kinetic energy (Ej(in) 
for the parameters characterizing the QR observed at in- 
termediate bias. As it is clear form FigfMl (Eja n ) can 
decrease significantly at intermediate Vbias- The effect 
becomes less and less evident decreasing the intermolec- 
ular hopping and disappears at t — 0.2. It is interesting 
to note that the potential energy curves show almost the 
same trend (see FigfPfl panel A). Therefore, for suffi- 
ciently large t and small E p , the oscillator overall energy 
decreases as a function of bias voltage. 

The behaviour of the average energy of the oscillator 
as function of bias voltage is determined by net balance 
of energy exchanged by the junction: after an increasing 
trend in the small bias regime, where the energy pumped 
by the bias exceeds that ceased to the electrons by the 
oscillator, the decreasing behaviour in the intermediate 
bias regime is due to the opposite physical mechanism: 
the energy ceased to the electron system by the oscillator 
exceeds that pumped by the bias. 

This v v transition ' ' occurs for that particular range 
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of bias voltages where the molecular energy levels are 
going through the bias window, with a resulting strong 
current enhancement (electronic resonance). In partic- 
ular, when the electron molecular levels enter the bias 
window completely, in the case of symmetric bias unbal- 
ance and for E g = 0, we expect that the electronic con- 
ductance reaches its maximum. Remarkably, comparing 
FigfT5land FigJT¥l one can observe that the conductance 
maxima correspond to kinetic energy minima, shifted by 
a quantity close to the EOC strength a. Physically, as 
a consequence of the SSH coupling with the oscillator, 
the current enhancement is followed by a strong effec- 
tive absorption of energy of the electron system from the 
oscillator. 



2. Electronic transport properties 

In order to evaluate the current through the molec- 
ular system in SSH model, we use the Meir-Wingreen 
formula^ for non interacting molecular levels, special- 
ized to our two-level case 

I{x) = | J - f R (u))Tr{G a r L G r T R }, 

(44) 

where the matrices / R are given by 

KT ( 1 l\ _ KT ( 1 -1\ 

v ^ = -\x ij' r * = T(-i i J' (45) 

and bold G r ' a indicate retarded (advanced) matrix green 
functions fEq|35l). We have explicitly indicated that the 
current depends on the deformation x of the oscillator, 
so that it has to be averaged over the probability distri- 
bution function P(x). 

Here, we focus on two particular physical regimes, pre- 
viously discussed in the analysis of the static approxi- 
mation: the weak coupling (E p « hT) and the strong 
coupling (KT « E p ) limits, varying arbitrarily the in- 
termolecular energy scale t. As we shall see, in both 
regimes, the direct coupling of the electron-oscillator in- 
teraction to the intermolecular hopping makes the role of 
the dynamical fluctuations crucial to determine correct 
results. 

In the weak coupling regime, it is interesting to ob- 
serve that, as in the AH model, the dynamical correc- 
tions renormalize and broaden the electronic resonances 
fFigfTS)) with respect to the static solution. In particu- 
lar, in panel A of FigfTSl we note that the static approx- 
imation exceeds the maximum value of non interacting 
conductance and shows a region of negative conductance 
at intermediate bias. However, when the dynamical con- 
tributions are included (square (green) curve), the effect 
on the conductance is dramatic washing out all the struc- 
tures observed in the static approximation. More inter- 
esting are the cases of panels B-C-D of FigfTSl where 
again the static approximation shows the spurious result 
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FIG. 15. (Color online). Panel A-D: Conductance as function 
of the bias at ujq/T = 0.1, E g = 0.0, interaction strength E p = 
0.2, for different values of intermolecular hopping t = 0.2 — 
0.6 — 1.0 — 2.0. The solid (black) line indicates non interacting 
curve. The dashed (red) line and the square (green) lines refer 
to the static and the dynamical approximation, respectively. 
The value of too shown in the figure is in F units, while all 
other quantities (E g , E p , t and eVuas) are expressed in hT 
units. 

of conductance greater than 2Go, while the dynamical ap- 
proximation renormalizes and broadens the peak of con- 
ductance to bias values where the static approximation 
shows small electric conduction. Even in the weak cou- 
pling regime, the inclusion of the dynamical fluctuations 
is essential to obtain correct results for the electronic con- 
duction. 

Finally, we examine the electronic transport proper- 
ties in the strong coupling regime {E p » hT), where 
moreover SP >> t. In this case, we expect strongly non 
linear behaviour of I-Vs in the infinite mass (static) limit 
for the oscillator. In Fig[T()]we show the current voltage 
characteristic for strong interaction, E p = 1.4, at fixed 
adiabatic ratio ojq/T = 0.1, gate voltage E g — 2.0 and for 
different small values of intermolecular hopping t < < hT 
(t = 0.15 (black) dashed, t = 0.2 (red) dashed dotted, 
t = 0.25 (blue) short dashed dotted line). In panel A 
we show a comparison between the non interacting and 
static approximation. The static approximation shows 
an interesting region of Negative Differential Resistance 
(NDR), as a consequence of the rich structure of the min- 
ima of the generalized potential described in the previous 
subsection (see also FigflTl Panel B). At intermediate 
bias voltage, a strong current currying region appears. 
This corresponds to 

) ~ ( n 72) — ^0.5 for which 
the electronic levels renormalize in the bias window with 
an effective energy larger than the " bare ' value. Then, 
for sufficiently large bias, the minimum corresponding to 
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FIG. 16. (Color online). Panel A: Current- Voltage charac- 
teristic in the static approximation for cjo/T = 0.1, E g = 2 
and strong coupling E v — 1.4 for different values of inter- 
molecular hopping (t = 0.15 dashed (black), t = 0.2 dashed 
dotted (red), t — 0.25 short dashed dotted (blue)). The non 
interacting quantities (t = 0.15 solid (black), t = 0.2 dotted 
(red), t = 0.25 short dotted (blue)) are also shown. Panel B: 
Current- Voltage characteristic in the dynamical approxima- 
tion for the same parameter values of panel A. The value of 
ojo shown in the figure is in V units, while all other quantities 
(Eg, E p , t and eVuas) are expressed in hY units. 



x ~ prevails, determining a strong current reduction 
due to the drop of the hoppings to their non interact- 
ing "bare' values. As in the case described above, one 
can note (panel B) that the dynamical corrections wash 
out all the features of the static approximation. There 
is a very small conduction threshold after which one not 
observe NDR features. Again, we observe that the in- 
clusion of dynamical corrections are very important for 
a correct description of the SSH model while the static 
approximation can easily lead to erroneous conclusions. 



IV. CONCLUSIONS AND DISCUSSIONS 

In this paper we have derived and studied the stochas- 
tic Langevin equation for the dynamics of an oscilla- 
tor mode coupled to a voltage-biased molecular junction 
in the adiabatic limit. Using the generalization of the 
Keldysh formalism to time dependent cases, we were able 
to show, in agreement with other approaches, that the os- 
cillator dynamics is controlled by an effective potential as 
well as by damping and fluctuating terms coming from 
the time depending electronic Green function. Actually, 
we have built up an expansion for the molecular level 
Green function in the velocity of the oscillator. In this 



way we have shown that the quantum effects, "hidden' 
in the stochastic equation, come only from the electronic 
subsystem. Solving numerically the Langevin equation, 
we have calculated the position and velocity distribution 
probabilities of the classical oscillator. We have focused 
our attention on the properties of the velocity distribu- 
tion function showing that, for sufficiently large bias volt- 
age, it loses its gaussian character. In addition, we have 
established the range of validity of the adiabatic approx- 
imation underlying the stochastic approach by setting 
QR (small bias), CAR (intermediate bias), CNAR (large 
bias) regimes. The criterion is based on the compari- 
son of the average kinetic energy of the oscillator with 
the Debye temperature (&_bTd ~ Uuio) to distinguish be- 
tween QR against CAR regimes, and with electron energy 
scale (~ hT) to distinguish between CAR against CNAR 
regimes. 

We applied our analysis to two simple models of 
molecules. 

For the single site AH model, the analysis of the va- 
lidity of the adiabatic approximation has allowed us to 
build up a diagram for the validity of classical approxima- 
tion in the plane (E p -Vbias), showing that the quantum 
effects are relevant only in a very narrow region if the adi- 
abatic ratio is smaller than all other energy scales. More- 
over, we have studied the current- voltage characteristic 
and the conductance, observing a dynamical reduction of 
the polaronic shift and the broadening of the electronic 
resonance due to the average on the nonequilibrium po- 
sition distribution probability of the oscillator. In the 
non-gaussian intermediate bias regime and for sufficiently 
large interaction strength, the kinetic energy shows an in- 
teresting non monotonic behaviour. Correspondingly, we 
observe in the transport properties a strong enhancement 
of conduction with respect to the infinite mass approxi- 
mation (static limit). 

We have also studied the case of a molecular Hamil- 
tonian composed by a couple of sites interacting with a 
single vibrational mode in the SSH model. In this case, 
because of the direct coupling of the electron-oscillator 
interaction to the intermolecular hopping, the role of the 
dynamical fluctuations becomes crucial to determine the 
physical scenario described by the model. The new in- 
termolecular electronic hopping energy scale t introduces 
a reduction of the CAR in the validity diagram. The 
new feature is the occurrence of small QRs for suffi- 
ciently small coupling E p , at intermediate bias voltages. 
For strong enough electron-phonon coupling E p , these re- 
gions disappear. In this region of parameters the average 
dynamical kinetic energy decreases as the bias voltage 
increases. Also the potential energy curves show this 
behaviour. Therefore, the oscillator overall energy de- 
creases as a function of bias voltage. This loss of energy 
occurs for that particular range of bias voltages where the 
molecular energy levels enter in the bias window. Corre- 
spondingly as in the AH model case, we observe in the 
transport properties an enhancement of conduction with 
respect to the infinite mass approximation. Remarkably, 
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in the case of small electron-oscillator interaction and count the coupling with the leads, 



for Eg = 0, we found that the maxima of conductance 
correspond to the minima of the kinetic energy, shifted 
by the EOC strength a. Finally, within this model, the 
dynamical corrections on the transport properties can- 
cel out completely the " detailed ' features (like NDR) 
present in the static case. As main result, we can con- 
clude observing that the inclusion of dynamical effects of 
the oscillator motion strongly modifies the physical sce- 
nario which would be obtained by a static description, 
even if the oscillator dynamics is much slower than the 
electron tunneling rate. 

We end this section noting that it could be of out- 
standing interest to study the possibility to include the 
quantum correction to the oscillator dynamics in the 
small bias regime classified as Quantum Region (QR). In 
this direction, Millis et a/. 40 find in the quasi-equilibrium 
regime E p >> luq >> Vbias a quantum contribution to 
the effective temperature of the oscillator in addition to 
the diffusive one. At finite mass m, nearby the " gaus- 
sian fluctuation ' paths involving small excursions (char- 
acteristic frequency ujq) from the minima of the static 
potential, quantum tunneling processes become impor- 
tant. The inclusion of the quantum corrections in our 
approach, within the minimal models considered, is un- 
der investigation. 
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Appendix A: Abiabatic Approximation 

Here, we show how the adiabatic approximation on the 
electronic Green function Eq.® works. 

In order to implement the adiabatic approximation, it 
is convenient to write the Dyson equation for the molec- 
ular retarded Green function (Eq.((9])) 



G r (t,t') = g r {t,t') + J dtxj dt 2 G r (t,t!) 

+ ^ el _ ph (t 1 ,t 2 )g r (t 2 ,t% (Al) 

where the Green function g r (t,t') already takes into ac- 



g^{t,t , ) = - l -e{t-t')e-<^- lV ^ t - t '\ (A2) 

Now, we reparametrize the retarded electron-oscillator 
self-energy separating slow and fast times scales (in the 
following, for sake of simplicity, we drop the label e i-ph 
of the self-energy) 



(A3) 



According to the approach used in Ref.26, we expand 

i+i 
2 



Eq. (|A3[) with respect to the slow mean time tl t t2 about 



a generic time to belonging to the interval [t, t'] 



E 1 



U -h 



(A4) 



with 



S5(t , h - t 2 ) = XxitoWh - t 2 ) (A5) 

S[ (to, - *a) = (^y^ - ta ) A *(*o)<% - h). 

(A6) 

The adiabatic expansion 

G r {t, t')~G r (t ,t-t') + Gl(t , t - t') (A7) 

for the Green function follows from that for the self- 
energy via the Dyson equation Eq. (|Al[) . We can now 
introduce the Fourier transforms G^^^o^uj) — J d(t — 
t ,yiu{t-t')/h G r a/ ^ t _ t /y Since Qur gQal ig an adia _ 

batic expansion of the electronic obscrvables at time t, we 
choose to = t. One can easily show that this is the only 
choice able to recover the fluctuation-dissipation theorem 
at vanishing bias voltage (equilibrium condition) for the 
Langevin equation (Eq ll7[) we derive in Section II.B.l. 
From the Dyson equation Eq. (|Aip . taking into account 
Eq.dA"5l) and Eq. (fA%)) . we thus find 



G r (Lu,t)~G r (u,t) + G r 1 (co,t) 



with 



G r (t,io) 



huj - E g (t) + ihr /2' 



(A8) 



(A9) 



GU^) = z^^^G S (^ ) , (A10) 

obtaining a correction which is linear in the velocity of 
the oscillator ^ = A§f . 
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